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! Abstract 

A non-vanishing Lyapunov exponent Ai provides the very definition of deter- 
ministic chaos in the solutions of a dynamical system, however no theoretical 

: 

mean of predicting its value exists. This paper copes with the problem of 

> : 

■ analytically computing the largest Lyapunov exponent Ai for many degrees 
C$ • of freedom Hamiltonian systems as a function of e = E/N, the energy per 

degree of freedom. The functional dependence Ai(e) is of great interest be- 
cause, among other reasons, it detects the existence of weakly and strongly 
chaotic regimes. This aim - analytic computation of Ai(e) - is successfully 
reached within a theoretical framework that makes use of a geometrization 
of newtonian dynamics in the language of Riemannian differential geometry. 
A new point of view about the origin of chaos in these systems is obtained 
independently of the standard explanation based on homoclinic intersections. 
Dynamical instability (chaos) is here related to curvature fluctuations of the 
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manifolds whose geodesies are natural motions and is described by means of 
Jacobi - Levi-Civita equation (JLCE) for geodesic spread. In this paper it is 
shown how to derive from JLCE an effective stability equation. Under general 
conditions, this effective equation formally describes a stochastic oscillator; an 
analytic formula for the instability growth-rate of its solutions is worked out 
and applied to the Fermi-Pasta-Ulam /3-model and to a chain of coupled rota- 
tors. An excellent agreement is found between the theoretical prediction and 

numeric values of Ai(e) for both models. 
PACS numbers(s): 05.45.+b; 02.40.-k; 05.20.-y 
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I. INTRODUCTION 



During the last two decades or so, there has been a growing evidence of the indepen- 
dence of the two properties of determinism and predictability of classical dynamics. In fact, 
predictability for arbitrary long times requires also the stability of the motions with respect 
to variations - however small - of the initial conditions. 

With the exception of integrable systems, the generic situation of classical dynamical 
systems describing, say, N particles interacting through physical potentials, is instability 
of the trajectories in the Lyapunov sense. Nowadays such an instability is called intrinsic 
stochasticity, or chaoticity, of the dynamics and is a consequence of nonlinearity of the 
equations of motion. 

Likewise any other kind of instability, dynamical instability brings about the exponential 
growth of an initial perturbation, in this case it is the distance between a reference trajectory 
and any other trajectory originating in its close vicinity that locally grows exponentially in 
time. Quantitatively, the degree of chaoticity of a dynamical system is characterized by 
the largest Lyapunov exponent Ai that - if positive - measures the mean instability rate of 
nearby trajectories averaged along a sufficiently long reference trajectory. The exponent Ai 
also measures the typical time scale of memory loss of the initial conditions. 

Let us remember that if 

x i = X\x l ...x N ) (1) 

is a given dynamical system, i.e. a realisation in local coordinates of a one-parameter group 
of diffeomorphisms of a manifold M, that is of 0* : M — > M, and if we denote by 

c = r k [x{t)\e (2) 

the usual tangent dynamics equation, i.e. the realisation of the mapping dcf) 1 : TM — > TM, 
where TM is the tangent bundle of M and is the Jacobian matrix of [X*], then the 
largest Lyapunov exponent Ai is defined by 
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and, by setting A[x(t),£{t)] = ? JW)}^ = ft/ft = kfM^l this can be for- 
mally expressed as a time average 

A 1 = lim 1 f drA[x(r),ar)} ■ (4) 

Even though Ai is the most important indicator of chaos of classical [ij dynamical systems, 
it is used only as a diagnostic tool in numerical simulations. With the exception of a few 
simple discrete-time systems (maps of the interval), no theoretical method exists to compute 
Ai @. This situation reveals that a satisfactory theory of deterministic chaos is still lacking, 
at least for systems of physical relevance. 

In the conventional theory of chaos, dynamical instability is caused by homoclinic in- 
tersections of perturbed separatrices, however this theory has many problems: i) it needs 
action-angle coordinates, ii) it works in conditions of weak perturbation of an integrable sys- 
tem, in) to compute quantities like Mel'nikov integrals one needs the analytic expressions 
of the unperturbed separatrices: at large iV this is hopeless, moreover the generalization 
of Poincare-Birkhoff theorem is still problematic at iV > 2; iv) finally, there is no compu- 
tational relationship between homoclinic intersections and Lyapunov exponents. Therefore 
this theory seems not adequate to treat chaos in Hamiltonian systems with many degrees of 
freedom at arbitrary degree of nonlinearity, with potentials that can be hardly transformed in 
action-angle coordinates, not to speak of accounting for phenomena like the transition from 
weak to strong chaos in Hamiltonian systems 0,|J. Motivated by the need of understanding 
this transition from weak to strong chaos, we have recently proposed |5|-[T0| to tackle Hamil- 



tonian chaos in a theoretical framework different from that of homoclinic intersections. This 
new method makes use of the well-known possibility of formulating Hamiltonian dynamics 
in the language of Riemannian geometry so that the stability or instability of a geodesic 
flow depends on curvature properties of some suitably defined manifold. 

In the early 1940s, N. S. Krylov already got a hold of the potential interest of this 
differential-geometric framework to account for dynamical instability and hence for phase 



space mixing |TTJ. The follow-up of his intuition can be found in abstract ergodic theory 
fl"2|j and in a very few mathematical works concerning the ergodicity of geodesic flows of 



physical interest [|T^,|T4|]. However, Krylov's work did not entail anything useful for a more 
general understanding of chaos in nonlinear newtonian dynamics because one soon hits 
against unsurmountable mathematical obstacles. By filling certain mathematical gaps with 
numerical investigations, these obstacles have been overcome and a rich scenario emerged 
about the relationship between stability and curvature 

Based on the so-obtained information, the present paper aims at bringing a substantial 
contribution to the development of a Riemannian theory of Hamiltonian chaos. The new 
contribution consists of a method to analytically compute the largest Lyapunov exponent 
Ai for physically meaningful Hamiltonian systems of arbitrary large number of degrees of 
freedom. A preliminary and limited account of the results presented here can be found in 
Ref. 0. 

The paper is organized as follows: Section [TT] is a sketchy presentation of the geometriza- 
tion of newtonian dynamics; Section |inj contains the derivation of an effective stability 
equation from Jacobi - Levi-Civita equation for geodesic spread and an analytic formula for 
Ai; Section [TV] contains the application of the general result to the practical computation of 
Ai in the Fermi-Pasta-Ulam /3-model and in a chain of coupled rotators. Some concluding 
remarks are presented in Section [V]. 



II. GEOMETRIZATION OF NEWTONIAN DYNAMICS 

Let us briefly recall how newtonian dynamics can be rephrased in the language of Rie- 
mannian geometry. We shall deal with standard autonomous systems, i.e. described by the 
Lagrangian function 

C — T — V — -a^j - Vfa, ...,q N ) , (5) 
so that the Hamiltonian function 7i = T + V = Eisa, constant of motion. 
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According to the principle of stationary action - in the form of Maupertuis - among all 
the possible isoenergetic paths j(t) with fixed end points, the paths that make vanish the 
first variation of the action functional 

A= Pidqi= -wrQidt (6) 

J<y(t) J-y(t) CQi 

are natural motions. 

As the kinetic energy T is a homogeneous function of degree two, we have IT = qidC/dqi , 
and Maupertuis' principle reads 



5A = 5 2Tdt = . (7) 

y 7 (t) 

The configuration space M of a system with N degrees of freedom is an N- dimensional 
differentiable manifold and the lagrangian coordinates (qi, ■ ■ ■ ,Qn) can be used as local 
coordinates on M. The manifold M is naturally given a proper Riemannian structure. In 
fact, let us consider the matrix 

9ij = 2[E - V{q)] aij (8) 

so that (0) becomes 

6 f 2Tdt = 5 [ (^<f<f) 1/2 dt = 5 [ ds =0 , (9) 

thus natural motions are geodesies of M, provided we define ds as its arclength. The metric 
tensor gj of M is then defined by 

9. J = 9ijdq i ®dq j (10) 

where (dq l , . . . , dq N ) is a natural base of T*M - the cotangent space at the point q - in the 
local chart (q 1 , . . . , q ). This is known as Jacobi (or kinetic energy) metric. Denoting by V 
the canonical Levi-Civita connection, the geodesic equation 

V^7 = (11) 
becomes, in the local chart (q 1 , . . . , q N ), 



^v +r i^!^! = 0) (12) 

ds 2 3 ds ds 

where the Christoffel coefficients are the components of V defined by 

r }fc = ( d 9*» V i efc ) = (%fcm + 9 k9mj ~ d m g jk ) , (13) 

where di = d/dq l . Without loss of generality consider = 2[E — V(q))6ij, from Eq. ([12]) 
we get 



d 2 q i 1 



+ 



' d(E-V)dq j dq i ij d{E-V) dq k dq m 
2 — ^ : — ; — g J — ^ gk 



, (14) 



ds 2 2(E — V) [ dqj ds ds dqj ds ds 

and, using ds 2 = 2(E — V) 2 dt 2 , we can easily verify that these equations yield 

d ^ dV 

-W = -dq- ^ = 1 '-' iV - ( 15 ) 

which are Newton equations. 

As already discussed elsewhere |5]||, there are other possibilities to associate a Rieman- 
nian manifold to a standard Hamiltonian system. Among the others we mention a structure, 



defined by Eisenhart |15|, that will be used in the following for computational reasons. In 
this case the ambient space is an enlarged configuration space-time Mxl 2 , with local coordi- 
nates (q°, q 1 , . . . , q N , q N+1 ), with (q 1 , . . . , q N ) G M, q° G E is the time coordinate, q N+1 G M 
is a coordinate closely related to the action; Eisenhart defines a pseudo-Riemannian non- 
degenerate metric g E on M x M 2 as 

ds 2 E = g^ dq* 1 <g> dq v = a {j dq 1 <g> dq j - 2V(q) dq° ® dg° + dg° ® rfg 7V+1 + rfg Af+1 ® . (16) 

Natural motions are now given by the canonical projection ir of the geodesies of (M x IR 2 , g^) 
on configuration space-time: 7r : M x M 2 — > M x M. However, among all the geodesies of gg 
we must consider only those for which the arclength is positive definite and given by 

ds 2 = g^dq^dq" = 2C 2 dt 2 , (17) 

or, equivalently, we have to consider only those geodesies such that the coordinate q N+1 
evolves according to 



g N + l = C 2 t + C 2_ r CdTj (1 8 ) 
JO 

where C and C\ are real constants. Since the values of these constants are arbitrary, we fix 
C 2 = 1/2 in order that ds 2 = dt 2 along a physical geodesic. For a diagonal kinetic energy 
matrix = 5ij, the non vanishing components of the connection V are simply 

Too = -r£ +1 = diV , (19) 

therefore it is easy to check that also the geodesies of g E yield Newton equations together 
with the differential versions of Eq. fll8D and of q° — t (details can be found in ||[| ) . 

III. GEOMETRIC DESCRIPTION OF DYNAMICAL INSTABILITY 

The actual interest of the Riemannian formulation of dynamics stems from the possibil- 
ity of studying the instability of natural motions through the instability of geodesies of a 
suitable manifold, a circumstance that has several advantages. First of all a powerful math- 
ematical tool exists to investigate the stability or instability of a geodesic flow: the Jacobi 
- Levi-Civita equation (JLC) for geodesic spread. The JLC equation describes covariantly 
how nearby geodesies locally scatter and it is a familiar object both in Riemannian geometry 
and theoretical physics (it is of fundamental interest in experimental General Relativity). 
Moreover the JLC equation relates the stability or instability of a geodesic flow with curva- 
ture properties of the ambient manifold, thus opening a wide and largely unexplored field of 
investigation of the connections among geometry, topology and geodesic instability, hence 
chaos. 

A. Jacobi - Levi Civita equation for geodesic spread 

A congruence of geodesies is defined as a family of geodesies {7 T (s) = 7(s,r) |r e 1} 
that, originating in some neighbourhood X of any given point of a manifold, are differentiably 
parametrized by some parameter r. Choose a reference geodesic 7(5, To), denote by j(s) the 
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field of vectors tangent at s to 7 and denote by J(s) the field of vectors tangent at r to the 
curves 7 s (r) at fixed s. The field J = (d^/dr)^ is known as geodetic separation field and 
it has the property: C^J = 0, where L is the Lie derivative. Locally we can measure the 
distance between two nearby geodesies by means of J. 

The evolution of the geodetic separation field J conveys information about stability or 
instability of the reference geodesic 7, in fact, if || J\\ exponentially grows with s then the 
geodesic is unstable in the sense of Lyapunov, otherwise it is stable. 

The evolution of J is described by fl9| 



^l + R^(s) } J(s))^s) = 0, (20) 

known as Jacobi - Levi-Civita (JLC) equation. Here J(s) G T 7 ( S )M; R(X, Y) = V^Vy — 
VyVx — V[x,y] is the Riemann-Christoffel curvature tensor; 7 = d^/ds; V/ds is the covari- 
ant derivative and 7(5) is a normal geodesic, i.e. such that s is the length. In the following 
we assume that J(s) is normal, i.e. (J, 7) = 0. This equation relates the stability or instabil- 
ity of nearby geodesies to the curvature properties of the ambient manifold. If the ambient 
manifold is endowed with a metric (e.g. Jacobi or Eisenhart) derived from the Lagrangian 
of a physical system, then stable or unstable (chaotic) motions will depend on the curvature 
properties of the manifold. Therefore it is reasonable to guess that some average global 
geometric property will provide information, at least, about an average degree of chaoticity 
of the dynamics independently of the knowledge of the trajectories, that is independently of 
the numerical integration of the equations of motion. 
In local coordinates the JLC equation (|20"D reads as 

^ + R ^ J (21) 

where R l j k i = (dq\ R(e^)i e (i)) e (j)) are the components of the curvature tensor, and the 
covariant derivative is (VJ l /ds) = dJ l /ds + T l j k J k dq^ /ds. There are 0(N 4 ) of such compo- 
nents, iV = dim M, therefore - even if this number can be considerably reduced by symmetry 
considerations - equation (|2lD appears untractable already at rather small N. It is worth 



mentioning that some exception exists. Such is the case of isotropic manifolds for which 



can be reduced to the simple form 

V 2 J { 



+ KJ' l = i = l,...,N, (22) 



ds 2 

where K is the constant value assumed throughout the manifold by the sectional curvature. 

The sectional curvature of a manifold is the iV-dimensional generalization of the gaussian 
curvature of two-dimensional surfaces of IR 3 . Consider two arbitrary vectors 1,7 e T X M, 
where x G M is an arbitrary point of M, and define 

\\X AY\\ = (||X|| 2 ||F|| 2 - {X,Y)f 2 (23) 

if ||X A Y || ^ the vectors X, Y span a two-dimensional plane 7r C T x M, then the sectional 
curvature at x relative to the plane ir is defined by 

KlXiY) = KM= wgM (24) 

which is only a property of M at x independently of X, Y G 7r (Gauss' theorema egregium). 
For an isotropic manifold K(x, ir) is also independent of the choice of ir and thus, according 
to Schur's theorem, K turns out also independent of x G M. 
Unstable solutions of the equation (|22]) are of the form 



J(s) = w(0)(-K)- 1/2 sinh (y^Ks^j , (25) 

once the initial conditions are assigned as J(0) = and dJ(0)/ds = w(0) and K < 0. In 
abstract ergodic theory geodesic flows on compact manifolds of constant negative curvature 
have been considered in classical works [TB| . In this case the quantity \J —K - uniform on 
the manifold - measures the degree of instability of nearby geodesies. 

While Eq. ( p2j) holds true only for constant curvature manifolds, a similar form of general 
validity can be obtained for JLC equation at N = 2. 

In this low-dimensional case Eq. (|2l| ) is exactly rewritten as 

§ + l -K(s)J = 0, (26) 
10 



where a parallely transported frame is used and 1Z(s) is the scalar curvature. Using Jacobi 
metric one finds (N = 2): 11 = AV/W 2 + (VV) 2 /W 3 , with W = E - V, so that for smooth 
and binding potentials TZ can be negative only where AV < 0, i.e nowhere for nonlinearly 
coupled oscillators as described by the Henon-Heiles model or for quartic oscillators ||10| . 
AV < is only possible if the potential V has inflection points. 

Recent detailed analyses of two-degrees of freedom systems |9|jl(| have shown that chaos 
can be produced by parametric instability due to a fluctuating positive curvature along the 
geodesies. 

Let us remember that parametric instability is a generic property of dynamical systems 
with parameters that are periodically or quasi-periodically varying in time, even if for each 



value of the varying parameter the system has stable solutions [fT^ . A harmonic oscillator 
with periodically modulated frequency, described by the Mathieu equation, is perhaps the 
prototype of such a parametric instability mechanism. 

Numerical simulations have shown that all the informations about order and chaos ob- 
tained by standard means (Lyapunov exponent and Poincare sections) are fully retrieved by 
using Eq. (|2B|). As in the case of tangent dynamics, Eq. fl2"E| ) has to be computed along a 
reference geodesic (trajectory). 

Let us now cope with the large N case. It is convenient to rewrite the JLC equation (|21~D 
in the following form 

+ [RicfrOO, -VOO) As) ~ Ric(7(s), -7(a)) 700] + W{<y{s), J(s)) 7 (s) = , 

(27) 



where W is the Weyl projective curvature tensor whose components W l jkl are given by |L8 



W% = K 3kl - -^—^iA - R jk St) , (28) 

and Ric is the Ricci curvature tensor of components Rij = R m im j- Weyl's projective tensor 
W (not to be confused with Weyl's conformal curvature tensor) measures the deviation 
from isotropy of a given manifold. For an isotropic manifold W^ kl = 0, and we recognize 
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in (|27D equation (|22"D, in fact in this case Rjiq^q 1 /(N — 1) is just the constant value of 
sectional curvature. Remind that the Ricci curvature at x G M is K R (X^) = RjiX^X 1 ^ = 
J2a=i -^(a)) where , . . . , Xr N \ form an orthonormal basis of T X M. Hence we 



understand that Eq. (|2T| ) retains the structure of Eq. fl22| ) up to its second term that 
now has the meaning of a mean sectional curvature averaged, at any given point, over 
the independent orientations of the planes spanned by Xt a \ and Xn,)', this mean sectional 
curvature is no longer constant along 7(3). The last term of ([27]) accounts for the local 
degree of anisotropy of the ambient manifold. 

Let us now consider the following decomposition for the Jacobi field J 

J(s) = Y,Ms)e (t) (s) (29) 

i 

where {e(i) . . . e^)} is an orthonormal system of parallely transported vectors. In this ref- 
erence frame it is 

S^E^w < 30 > 

and the last term of ( f?7D is 

W^,J)^ = J2(W(i,J)i,e U) )e U) 
j 

= X}(W(7, Yl J * e «)7, e 0')> e U) ( 31 ) 
j * 

= ^(H/(7,e (i) )7,e (i) ) Jie (i) , 
y 

the same decomposition applies to the third term of Eq. ( p7|) which is finally rewritten as 

+ fcfl ( s ) j. + ^ ( w . . + r . .) Jt = (32) 



where fcjj = K R /(N - 1), u^- = (^(7, e (i) )7, e (i) ) and r y = (Ricft, e (i) )7, e^)/(N - 1). 
Of course /cr is independent of the coordinate system. The elements to^ still depend on 
the dynamics and on the behavior of the vectors eo.)(s), thus, in order to obtain a stability 
equation, for the geodesic flow, that depends only on average curvature properties of the 
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ambient manifold, we try to conveniently approximate the Wij. To this purpose define at 
any point x G M the trilinear mapping R' : 7 ,. \ I x T X M x T X M —>■ T X M by 

(R'{X, Y, U),Z) = (X, U) (Y, Z) - (Y, U) (X, Z) (33) 



for all X, Y, U, Z G T X M. It is well known [Tj| that, if and only if M is isotropic then 
R = KqR', where R is the Riemann curvature tensor of M and Kq is the constant sectional 
curvature. 

Let us now assume that the ambient manifold is quasi-isotropic, i.e. that it looks like an 
isotropic manifold after a coarse-graining that smears out all the metric fluctuations, and let 
us formulate this assumption by putting R ps K{s)R' and Ric « K(s)g, although K(s) is no 
longer a constant. Now we use (|33|) to find Wy ~ 8K(s)[(j,j)(e(i^,e(j^) — (e^,j)(j, e^)], 
then we use Ric oc g and g(j, J) — to find = thus Eq. ( p^ ) becomes 



+ fc fl (s) J,- + (Jif (s) Jj = , (34) 



by SK(s) = K(s) — K we denote the local deviation of sectional curvature from its coarse- 
grained value K, thus 5K(s) measures the fluctuation of sectional curvature along a geodesic 
due to the local deviation from isotropy. The problem is that 5K(s) still depends on a 
moving plane tt(s) determined by j(s) and J(s). In order to get rid of this dependence, 
consider that if x G M is an isotropic point then the components of the Ricci tensor are 
Rih = (N — l)K(x)gih and the scalar curvature is R = N(N — l)K(x); with these quantities 
one constructs the Einstein tensor Gih = Rih — \gihR whose divergence vanishes identically 
{Gih\i — 0) so that it is immediately found that, if a manifold consists entirely of isotropic 
points, then dK(x)/dx l = and so 8Kr(x) / 'dx l = 0, i.e. the manifold is a space of constant 
curvature (Schur's theorem |19|). Conversely, the local variation of Ricci curvature detects 
the local loss of isotropy, thus a reasonable approximation of the average variation SK(s) 
along a geodesic may be given by the variation of Ricci curvature. 

Next let us model 8K(s) along a geodesic by a stochastic process. In fact K(s) is ob- 
tained by summing a large number of terms, each one depending on different combinations 
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of the components of J and on the the coordinates q l , moreover, unless we tackle an inte- 
grable model, the dynamics is always chaotic and the functions q l (s) behave irregularly. By 
invoking a central-limit-theorem argument, at large N, SK(s) is expected to behave, in first 
approximation, as a gaussian stochastic process. More generally, the probability distribution 
V(SK) may be other than gaussian and in practice it could be determined by computing its 
cumulants along a geodesic 7(5). 

Now we make quantitative the previous statement - about using the variation of Ricci 
curvature along a geodesic to estimate SK(s) - by putting 

V(8K) ~ V(8K R ) . (35) 

Both 5K and 5Kr are zero mean variations, so the first moments vanish; according to (|35|) 
the following relation for the second moments will hold 

([K(s) -K] 2 ) s ~ j^([Kr(s) - <i^) s ] 2 >, , (36) 

where (-) s stands for proper-time average along a geodesic 7(5). Let us comment about the 
numerical factor in the r.h.s. of (|36| ) where a factor ™ might be expected. At increasing iV 
the mean square fluctuations of k,R drop to zero as because is the mean of independent 
quantities, however this cannot be the case of the mean square fluctuations of K, in fact 
out of the sum Kr of all the sectional curvatures, in Eq. ( [34]) only one sectional curvature 
is "picked-up" from point to point by 5K so that 5K remains finite with increasing N. 
Therefore, as the second cumulant of SK does not vanish with N, we have to keep finite 
the second cumulant of 8Kr, what is simply achieved by properly adjusting the numerical 
factor in Eq. fl36|). 

The lowest order approximation of a cumulant expansion of the stochastic process 5K(s) 
is the gaussian approximation 

5K(s) ~ -^=(5^)1/^(8) , (37) 

where r](s) is a random gaussian process with zero mean and unit variance. Finally, in order 
to decouple the stability equation from the dynamics, we replace time averages with static 
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averages computed with a suitable ergodic invariant measure fi. As we deal with autonomous 
Hamiltonian systems, a natural choice is the microcanonical measure on the constant energy 
surface of phase space PU| 

fi(x5(H- E) (38) 

so that Eq. ( |3~?D becomes 

6K(s)^-^=(5 2 K R yj* v (s) . (39) 



Similarly, k R (s) in Eq. fl3j|) is replaced by (k R )^, in fact at large N the fluctuations of k R - 
as already noticed above - vanish as i because the coarse-grained manifold is isotropic, so 
that we finally have 

|| + (A*)„V + -^=(S 2 K R )l/ 2 V (s)ij = , (40) 

where ip stands for any of the components J J , since all of them now obey the same effective 
equation of motion. The instability growth-rate of ip measures the instability growth-rate of 
|| J|| 2 and thus provides the dynamical instability exponent in our Riemannian framework. 
Equation © is a scalar equation which, independently of the knowledge of dynamics, pro- 
vides a measure of the average degree of instability of the dynamics itself through the behav- 
ior of ip(s). The peculiar properties of a given Hamiltonian system enter Eq. ([1(]) through 
the global geometric properties {k R )^ and (5 2 K R ) fl of the ambient Riemannian manifold 
(whose geodesies are natural motions) and are sufficient to determine the average degree of 
chaoticity of the dynamics. Moreover, according to (p8|), (k R )^ and (5 2 K R ) fl are functions 
of the energy E of the system - or of the energy density e = E/N which is the relevant 
parameter as N —>■ oo - so that from (^) we can obtain the energy dependence of the 
geometric instability exponent. 

B. An analytic formula for the largest Lyapunov exponent 



By transforming Eq. ([2t]) into Eq. ( |4"0"D the original complexity of the JLC equation 
has been considerably reduced: from a tensor equation we have worked out an effective 
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scalar equation formally representing a stochastic oscillator. In fact (|40|), with a self-evident 
notation, is in the form 

^ + Q(s)^ = (41) 

where Q(s) is a gaussian stochastic process. 

Now, passing from proper time s to physical time t, Eq. fl4l|) simply reads 

(Pip 

~d¥ 

where 



+ Q{t)ifj = , (42) 



m = (kn}, + -^=(5 2 K R ) l J* V (t) (43) 



if the Eisenhart metric is used [because of the affine parametrization of the arclength with 
time, Eq. flTTDI; if Jacobi metric is used, we have 



, WH ^ + (_M-) 2 + -(-)\ + i (44) 



[see Eq. (64) of [||] and Eq. (27) of ||], note that d/dt = q j (d/dq j ). Being interested in the 
large N limit, we replaced N — 1 with iV in Eqs. ( [43]) and ([4*4]) . Of course Ricci curvature 
has different expressions according to the metric used. 

The stochastic process Q(t) is not completely determined unless its time correlation 
function Tn(ti,t 2 ) is given. We consider a stationary and 5-correlated process Q(t) so that 
^n(ti,t 2 ) = r n (|t 2 - and 

T n (t)=ra 2 n 5(t) , (45) 

where r is a characteristic time scale of the process. In order to estimate r, let us notice 
that for a geodesic flow on a smooth manifold the assumption of 5-correlation of Q(t) will 
be reasonable only down to some time scale below which the different iable character of the 
geodesies will be felt. In other words, we have to think that in reality the power spectrum of 
Q(t) is flat up to some high frequency cutoff, let us denote it by v±\ therefore, by representing 
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the 5 function as the limit for v — ► oo of 8 v (t) = Bin j\ , a more realistic representation of 



the autocorrelation function in Eq. fl^o] ) could be r^(t) = cr n~ sm jj f * t ' > = r *o"n^*(*)) 

whence r* = l/V*. Notice that J^T^(t)dt = ra^ and f™Tn(t)dt = ^t*<Tq thus r = r*/2. 
For practical computational reasons it is convenient to use T^t) in the form given by Eq. 
( filf) (with the implicit assumption that v+ is sufficiently large), however, being finite, the 
definition r = will be kept. To estimate we proceed as follows. A first time scale, 
which we will refer to as r 1; is associated to the time needed to cover the average distance 



between two successive conjugate points along a geodesic pl| . In fact, at distances smaller 
than this one the geodesies are minimal and far from looking like random walks, whereas 
at each crossing of a conjugate point the separation vector field increases as if the geodesies 
in the local congruence were kicked (this is what happens when parametric instability is 



active). From Rauch's comparison theorem |19[ we know that if sectional curvature K is 



bounded as follows: < L < K < H , then the distance d between two successive conjugate 
points is bounded by ^7= < d < -^=. We need the lower bound estimate that, for strongly 
convex domains |[22| , is slightly modified to d > 
Hence we define T\ through 

\ds/ \ds/ 2y/Sh + <T 1 ; 

where (4^) is the average of the ratio between proper and physical time ((^) = 1 if Eisenhart 
metric is used) and the upper bound H of K is replaced by the iV-th fraction of a typical 
peak value of Ricci curvature, which is in turn estimated as its average Q plus the typical 
value 5K of the (positive) fluctuation, i.e. in a gaussian approximation 5K = a n . This 
time scale is expected to be the most relevant only as long as curvature is positive and the 
fluctuations, compared to the average, are small. 

Another time scale, referred to as r 2 , is related to local curvature fluctuations. These 
will be felt on a length scale of the order of, at least, / = l/y / o^' (the average fluctuation of 
curvature radius). The scale / is expected to be relevant one when the fluctuations are of the 
same order of magnitude as the average curvature. When the sectional curvature is positive 
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(resp. negative), lengths and time intervals - on a scale I - are enlarged (resp. shortened) 
by a factor (l 2 K/6) p5| , so that the period -J== has a fluctuation amplitude d 2 given by 

d 2 



i- /v .-'tt . re p] ac i n g ^ 1-,-y j^g m ost probable value one gets 



n - 



T2 



dt \ d _ / dt \ /2fi o 2n 
ds J 2 \ds J 6 v^o 



dt\ n 



1/2 




ds / Or 



(47) 



Finally r in Eq. ( [45]) is obtained by combining t\ with T2 as follows 

r- 1 = 2T- 1 = 2 (rf 1 + r^ 1 ) . (48) 
The present estimate of r is very close - though not equal - to the one of Ref . M . 



Whenever Q(t) in Eq. (J42I) has a non- vanishing stochastic component the solution ip(t) 
has an exponentially growing envelope [^4| whose growth-rate provides a measure of the 
degree of chaoticity. Let us call this quantity Lyapunov exponent and denote it by A. In the 
next Section we shall make more precise the relationship of A with the conventional largest 
Lyapunov exponent. 

Our exponent A is defined as 



x=lim L log m±m. 

t^2t 6 ^ 2 (0)+^ 2 (0) 



(49) 



where if>(t) is solution of Eq. (f£^). 

The ratio (ilj 2 (t)+ip 2 (t))/(ip 2 (0)+ip 2 (0)) is computed by means of a technique, developed 
by Van Kampen and sketched in Appendix A, which is based on the possibility of computing 
analytically the evolution of the second moments of if) and if), averaged over the realizations 
of the stochastic process, from 



d 
dt 





( 








2 


\ 


' (ip 2 ) N 


m - 

















\ 




1 





/ 





(50) 



where f2 and <7q are respectively the mean and the variance of Q(t) above defined. By 
diagonalizing the matrix in the r.h.s. of floTf ) one finds two complex conjugate eigenvalues, 
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and one real eigenvalue related to the evolution of | y(4> 2 ) + (V> 2 )J- According to P§|) the 
exponent A is half the real eigenvalue. Simple algebra leads to the final expression 

A(O ,a n ,T) = ^(A--^y (51a) 



2 V 3A 



A= ^r+^/(^) 3 + (2^j . (51b) 

All the quantities Qq, <tq and r can be computed as static averages, therefore - within the 
validity limits of the assumptions made above - Eqs. Q5T|) provide an analytic formula to 
compute the largest Lyapunov exponent independently of the numerical integration of the 
dynamics and of the tangent dynamics. 

1. Lyapunov exponent and Eisenhart metric 

Let us consider dynamical systems described by the Lagrangian function (|5]) with a 
diagonal kinetic energy matrix, i.e. = 5ij, and let us choose as ambient manifold the 
enlarged configuration space-time equipped with the Eisenhart metric (Jl6|). 

Trivial algebra gives Tq = (dV/dqi) and = (—dV/dq 1 ) as the only nonvanish- 

ing Christoffel coefficients and hence the Riemann curvature tensor has only the following 
nonvanishing components 

Rm > = ■ (52) 



The JLC equation (|20|) is thus rewritten in local coordinates as 



ds ds 10-7 ds ds ds ds 

W p , (dq°\ 2 p t dq° dq> , dq> dq° 

dsds + ^ \ds + 00j ds ds + J0 ° ds ds 



V V TN+ 1 -i- n N + ldqi T° dqj -i- pN+i dqi p dq ° - n 

TsTs 3 + R ^ rf7 J Its' + ^° lb J ~ds~ - ° • (53) 



As = implies VJ°/ds = dJ°/ds and as R^ k = 0, we find that the first of these 
equations reads 
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d 2 7° 

hence J° does not accelerate and, without loss of generality, we can set J°(0) = J°(0) = 0, 
this yields (using VJ l /ds = df/ds + r ok q° J k + r ko q k J°) 

V 2 J* cPJ* 

~d^ = lt^ (55) 
and the second equation in (|53|) gives, for the projection in configuration space of the sepa- 
ration vector, 

d 2 J* d 2 V fdq°\ 2 

the third of equations fl55| ) describes the passive evolution of J^ 4 " 1 which does not contribute 
the norm of J because S'Tv+iAf+i = 0, so we can disregard it. 

As already mentioned in the previous Section, along the physical geodesies of g& it is 
ds 2 = (dq ) 2 = dt 2 therefore Eq. Q56D is exactly the usual tangent dynamics equation 
reported in the Introduction, provided that the obvious identification £ = (£ ? ,£ p ) = (J, J) is 
made. This clarifies the relationship between the geometric description of the instability of a 
geodesic flow and the conventional description of dynamical instability. It has been recently 
shown P,|ToH that the solutions of the equations (|56D and (p6|) (where 1Z is computed with 
Jacobi metric) are strikingly close one another in the case of two degrees of freedom systems. 
This result is reasonable because the geodesies of (M xl 2 ,^)- that are natural motions - 
project themselves onto the geodesies of (M,gj), and as the extra coordinates q° and q N+l 
do not contribute to the instability of the geodesic flow , both local and global instability 
properties must be the same with either Jacobi or Eisenhart metrics, independently of N. 

With Eisenhart metric the only nonvanishing component of the Ricci tensor is i?oo — 
AV, where A is the euclidean Laplacian in configuration space. Hence Ricci curvature is 
^r(q) = AV/(N — 1) (remember that we choose the constant C such that ds 2 = dt 2 along 



a physical geodesic) and the stochastic process f2(i) in ( f42|) is specified by 

n = (k R )^ = ^(AV) IM , (57a) 
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N 



(57b) 



2r 



(57c) 



2. Averages of geometric quantities 



Let us now sketch how to compute the mean and the variance of any observable function 
f(q), a geometric quantity of the chosen ambient manifold, by means of the microcanonical 
measure (]58"|), i.e. 

(/(?)>m = — / f(q)5(H(q,P)-E)dqdp (58) 
J 

where 



lu e = j 5(H(q,p) - E)dqdp (59) 

and q = (qi ■ ■ ■ (?jv)> P = (Pi ■ ■ - Pn)- By using the configurational partition function Zc(fl), 
given by 



Z c {(3)= / dqe 



-(3V(q) 



(60) 



where dq = Yli=i dq%, we can compute the Gibbsian average (f) G of the observable / as 



(f) G = [Zc(P)}- 1 / dqfMe-M*) . 



Whenever this average is known, we can obtain the microcanonical average of / 
following parametric form 



(61) 



in the 



(f)M = (f) G (P) 



(62) 
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By replacing / with the explicit expression for Ricci curvature kn = jjKr we can work 
out f^o- Notice that Eq. ( |62] ) is strictly valid in the thermodynamic limit; at finite A" it is 

At variance with the computation of (/), which is insensitive to the choice of the 
probability measure in the A" — > oo limit, computing the fluctuations of /, i.e. of 
{S 2 f) = jr((f — (f)) 2 ), by means of the canonical or microcanonical measures yields dif- 
ferent results. The relationship between the canonical - i.e. computed with the Gibbsian 
weight e~^ w - and the microcanonical fluctuations is given by the well known formula [27 

2 



f?_ \ d(f) G (P) 

Cy [ dp 



(63) 



where 



a 



V 



(5 2 d(E) 



(64) 



A^ dp 

is the specific heat at constant volume and P = P(e) is given in implicit form by the second 
equation in d62|). 



By replacing / with we can work out a^. 



IV. APPLICATIONS 

In this Section the Riemannian approach to Hamiltonian chaos described above is practi- 
cally used to compute Ai(e) for two different models: the Fermi-Pasta-Ulam (FPU) /5-model 
and a chain of coupled rotators. 

The choice of these models is motivated by the possibility of analytically computing, in 
the A^ — ► oo limit, the geometric quantities needed, and by their interest as mentioned in 
the following subsections. 



A. The Fermi-Pasta-Ulam /3-model 



The FPU /3-model is defined by the Hamiltonian |25 
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N 1 N 



(65) 



i=l i=l 

This is a paradigmatic model of nonlinear classical many-body systems that has been ex- 
tensively studied over the last decades and that stimulated remarkable developments in 
nonlinear dynamics, one example: the discovery of solitons. For a recent review we refer 
to [36|. Also the transition between weak and strong chaos has been first discovered in 
this model and then, the effort of understanding the origin of such a threshold has 
stimulated the development of the geometric theory presented here. 

Let us now compute the average Ricci curvature Sl and its fluctuations cx n . We have 
seen above that, using Eisenhart metric, k R is given by 

d 2 V(q) 



, if* 2 



1 dq. 



2 ' 



(66) 



for the FPU /3-model this reads 



N 



i - <k ) 



(67) 



8=1 



note that k R is always positive. 

In order to compute the Gibbsian average of k R and its fluctuations, we rewrite the 
configurational partition function as 

-oo N f N 

TT da,-, exn 1—3 \~~ ^ ,, , . 



Z c (a) 



] J dq,i exp 
i=i ^ i=i 

which, in terms of the arbitrary parameter a and of Z c , is expressed as Zc{ot 
Zc (a/3,[i/a) and leads to the following identity 

12/i 1 



(68) 



(k R )((3) = 2- 



I3N Z t 



c 



(69) 



a=l 



Thus we have to compute 

1 



NZ, 



c 





1 






Q =i"^ 



log Z c (a) 



(70) 



using 
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Z c (a) = [z c (a)] N f(a) , 



(71) 



where f(a) is a quantity 0(1), z c {a) is the single particle partition function |28| 



1\ ffifi 



-1/4 



z c (a) = r ( § ) V~f) exp(-a 2 e 2 )D_ 1/2 (at 



T is the Euler function, -D_i/2 is a parabolic cylinder function and 



(72) 



2fi 



1/2 



(73) 



The final result in parametric form of the average Ricci curvature of (M x M 2 , p s ) - with 
the constant energy constraint - is (details can be found in ||) 

3 D_3/ 2 (0) 



fto(e) -> < 



<fc*>(0) = 2 + 



8cr 



9 £Li /2 (0) 



3 1 g-gMg) 
2 0£>_ 1/2 (0)J 



(74) 



Let us now compute 



(75) 



iV N " " Jl/ ^~' AT' 

According to Eq. (0), first the Gibbsian average of this quantity, (5 2 kn) G (l3) = 
jj((Kr — (Kn)) 2 ) G (j3), has to be computed and then the correction term must be added. 



Now define 



N 



(76) 



after Eq. 



i-(S 2 K R ) G (P) = i-((K R -(K B )Y) 



2\G 



N 



N 



36/x 2 



((Q-(Q)n 



2\G 



(77) 



hence using Eq. (p 



logZ c (a) 



(78) 



and finally 
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N 



(S 2 K R ) 



G 



144/i 2 



da'- 



logz c (a) 



(79) 



a=l 



Simple algebra gives 
d 2 



da 2 



log z c (a) 



a=l 



4 I D„ 1/2 (9) 



D -3/2(0) 



(80) 



so that from Eq. ( [79]) we obtain 



(81) 



According to the prescription of Eq. fl63|), the final result for the fluctuations of Ricci 
curvature is 



e(9) 



N 



8// 



c v (6) V d/3 



(82) 



3 1 £>- 3/2 (e) 

9 2 e D. 1/2 (6) 



where (5 2 K R ) G (9) is given by fl8T|), the derivative part of the correction term is 

d(k R )(6) _ 3 9D 2 _ V2 (9) + 2(9 2 -l)D_ 1/2 (9)D_ V2 (9)-29D 2 _ 1/2 (9) 
d(3 8/i# 3 ^ 2 i /2 (^) 

and the specific heat per particle c y is found to be 

c v (0) = 16D 2 {e) K 12 + 202 ) D - x/M + 29D_ 1/2 (9)D_, /2 (9) 



(83) 



- 9 2 D_ 3/2 (9) [29D_ 1/2 (9)+D_ 3/2 (9)]} 



(84) 



The micro canonical averages in Eqs.(|74]) and fl32|) are compared in Figs. [I] and |2] with their 



corresponding time averages computed along numerical trajectories of the model (|65 ) at 
iV = 128 and N = 512 with \x = 0.1. The equations of motion are integrated using a third 
order bilateral symplectic algorithm which is a high precision numerical scheme. 

Though microcanonical averages are computed in the thermodynamic limit, the agree- 
ment between time and ensemble averages is excellent already at N = 128. 
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1. Analytic result for Ai(e) and its comparison with numeric results 



Now we use ([73]) and (|82|) to compute r according to its definition in (57c), then we 
substitute flo(e), cr^(e) and r(e) into Eq. (pT|) to obtain the analytic prediction for Ai(e) in 
the limit N — > oo. In Fig. ^ this analytic result is compared to the numeric values of Ai 
computed by means of the standard algorithm JJIJ at N — 256 and N = 2000 with fi = 0.1 
and at different e. The agreement between analytic and numeric results is strikingly good. 



B. A chain of coupled rotators 

Let us now consider the system described by the Hamiltonian 

N ( 2 ^ 

nip, q) = J2 \y + J[1 - cosfe+i - %)] r • (85) 

If the canonical coordinates ft and pi are given the meaning of angular coordinates and 
momenta, this Hamiltonian describes a linear chain of N rotators constrained to rotate on 
a plane and coupled by a nearest-neighbor interaction. 

This model can be formally obtained by restricting to one spatial dimension the clas- 
sical Heisenberg model whose potential energy is V — — JY2{ij) S« ' Sj> where the sum is 
extended only over nearest-neighbor pairs, J is the coupling constant and each Si has unit 
module and rotates on a plane. To each "spin" Si = (cos ft, sin q{) the velocity J^Sj = 
( — ^dt sm ^dt cos * s associated so that (|85|) follows from H = X^=i §S? — J j) Si • Sj. 

The Hamiltonian fl85|) has two integrable limits. In the limit of vanishing energy it 
represents a chain of harmonic oscillators 

N f 2 1 
H(p,q)^J2 | + AQi+i-Qi) 2 } > (86) 

whereas in the limit of indefinitely growing energy a system of freely rotating objects is 
found because of potential boundedness. 

The expression of Ricci curvature Kpt, computed with Eisenhart metric, is 
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i=l ft 1=1 

Let us observe that for this model a relation exists between potential energy V and Ricci 
curvature Kr. 

V(q) = JN-^. (88) 

This relation binds the fluctuating quantity that enters the analytic formula for Ai. This 
constraint does not exist for the sectional curvature thus a-priori it may be expected that 
some problem will arise. 

The configurational partition function for a chain of coupled rotators is 

n n r n 

Z C (P) = / I! dqi exp 1 ~P Yl J t 1 ~ cos (ft+i ~ ft)] 
J -"i=i { i=i 

N N 

= exp(-pJN) / Y[duJiexp(P J ^ cos uji) (89) 

J —IT ■ i ■ i 



1=1 1=1 

N (n„\N „(—\ 



= exp(-/5JiV)[/ ( / 3J)] w (27r) w ^(cJ) . 

where ioC^) — ^ Jq W e xcos9 d9 is the modified Bessel function of index zero; Ui = q i+1 — q iy 
i e (1, . . . , iV — 1),^ = ^ — q^, q = u depend on the initial conditions. The function g(cu) 
contributes with a term of O(jj) thus vanishing in the thermodynamic limit. 

In order to compute f2 and o"^ we follow the same procedure adopted for the FPU model, 
i.e. we define 

/•+*■ N f N } 

z c(a) = j J dq t exp <-(3^[l - a cos(g m - $)] ^ 

J -* i=i I i=i J 

= exp(-/3JJV) [7 (/3Jo ; )] JV ^(a;)(27r) iV (90) 

and by observing that 



AT/5 

we find Qq(e) m parametric form 



|;logZ c (a) 



(91) 



o=l 
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UP J) 
HP J) 



£ (/3) 



2/3 + v w);' 



(92) 



In order to work out the average of the square fluctuations of Ricci curvature we use the 
following identity 



da'- 



logZ c (a) 



(93) 



a=l 



whence 



i 2 G l2 (3JWJ) - hWJWM - pjiKpj) 
n {SKr) =4J WWJ) 



(94) 



The computation of the correction term 



d{k R )(fi) 
9P 



-\ 2 



I de ^ involves the following derivatives 



hp) 

dp 



i HP J) 
pJh(PJ) 



HP J) 
Hpj) 



(95) 



d(k R )(P) 



2J 2 



i HP J) 



HP J) 
HP J) 



op r pJHPJ) 

Finally, gluing together the different terms, we obtain 

4j pJiliPJ) - HPJ)HPJ) - PJitiPJ) 



(96) 



<P)=YP +J 



P il(pj) [i + 2 (pj) 2 ] - 2Pjhpj)hpj) - 2 [pjhpj)Y 

\ hpj)' 



HP J). ' 



(97) 



In Figs. |] and |5| the comparison between analytic and numeric results is provided for the 
average Ricci curvature and its fluctuations. The agreement between ensemble and time 
averages is very good. Time averages are computed along numerical trajectories of the 
model Hamiltonian ( |8~5D at iV = 150 and J = 1. The already mentioned high precision 
symplectic algorithm has been used also in this case. 
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1. Analytic result for Ai(e) and its comparison with numeric results 



By inserting into Eq. Q5T| ) the analytic expressions of VLq(e) and cr^(£) given in Eqs. 
and fl97p - and also r(e) which is a function of the latter quantities - we find Ai(e). 

In Fig. |6] the comparison is given between the analytic result so obtained and the outcome 
of numeric computations performed with the standard algorithm [ 30 1 . Figure |6] shows that 



there is agreement between analytic and numeric values of the largest Lyapunov exponent 
only at low and high energy densities. Likewise the FPU case, at low energy, in the quasi- 
harmonic limit, we find \\{s) oc e 2 . Whereas at high energy Xi(e) oc e" 1 / 6 , here Xi(e) 
is a decreasing function because at e — > oo the systems is integrable. In an intermediate 
energy range our theoretical prediction underestimates the actual degree of chaoticity of 
the system. It is worth mentioning that this energy range coincides with a region of fully 
developed (strong) chaos detected in this model by a completely different approach in Ref. 
H. In this already mentioned above - there was a-priori a reason to expect an 

inadequacy of the analytic prediction in some energy range. In fact, using Eisenhart metric, 
the explicit expression of the sectional curvature K(v,£) - relative to the plane spanned 
by the velocity vector v and a generic vector £_Lv (here we use £ to denote the geodesic 
separation vector in order to avoid confusion with J which is the notation for the coupling 
constant) - is 



dq° £* dq° £ fe d 2 V £*£ fe 



hence we get 



T * 

K{v, = cos(g m - «) [C +1 - C] 2 (99) 

for the coupled rotators model. We realize, by simple inspection of Eq. (|99|), that K can 
take negative values with non-vanishing probability regardless of the value of e, whereas - 
as long as e < J - this possibility is lost in the replacement of K by Ricci curvature that we 
adopted in our theory. In fact, because of the constraint flg8|), at each point of the manifold 
it is 
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k R (e) > 2(J - e) (100) 

thus our approximation fails in accounting for the presence of negative sectional curvatures 
at small values of e. In Eq. fl99|) the cosines have different and variable weights, — £ 1 ] 2 , 
that in principle make possible to find somewhere along a geodesic K < also with only 
one negative cosine. This is not the case of k R where all the cosines have the same weight. 
Therefore the probability of finding K < along a geodesic must be related to the probability 
of finding an angular difference greater than | between two nearest-neighboring rotators. If 
the energy is sufficiently low this event will be very unlikely, but we can guess that it will 
become considerable where the theoretical prediction is not satisfactory, i.e. when chaos is 
strong. Notice that the frequent occurrence of K < along a geodesic adds to parametric 
instability another instability mechanism that enforces chaos [Eq. (p5|)j . 

Our strategy is to modify the model for K(s) in some effective way that takes into 
account the mentioned difficulty of k^(s) to adequately model K(s). This will be achieved 
by suitably "renormalizing" Q or o~ n to obtain an effective gaussian process for the behavior 
of the sectional curvature. 

From Eq. ( |55| ) we see that N directions of the vector £ exist such that the sectional 
curvatures - relative to the N planes spanned by these vectors together with v - are just 
cos(gj + i — qi). Hence the probability P{e) of occurrence of a negative value of the cosine 
is used to estimate the probability of occurrence of negative sectional curvatures along the 
geodesies. This probability function has the following simple expression 

f n e(-cosx)e^ cos ^x ffe^^dx 

P(e) = — — = — (101) 

K) f^ePJcnvfa 2irI ({3J) ' 1 } 

where Q(x) is the Heavyside unit step function. 

The function P(s), reported in Fig. [7|, begins to increase at e ~ 0.2, just where the 
analytic prediction in Fig. || begins to fail, and when it approaches its asymptotic value of 
|, around the end of the knee, a good agreement is again found between theory and numeric 
results. The simplest way to account for the existence of negative sectional curvatures is to 
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shift the peak of the distribution V(5Kr) toward the negative axis. This is achieved by the 
replacement 

<M,)> d02) 

This correction neither has influence when P(e) ~ (below e ~ 0.2) nor when P(e) ~ 1/2 
(because in this case (k R (e)) — > 0). The value of the parameter a in (|102|) must be estimated 



a posteriori in order to obtain the best agreement between numerical and theoretical data 
over the whole range of energies. The result shown in Fig. || is obtained with a = 150, 
anyhow no particularly fine tuning is necessary to obtain a very good agreement between 
theory and numerical experiment. 

V. CONCLUDING REMARKS 

The present paper contains a substantial progress along the research line initiated in Ref. 
H where it was proposed to tackle Hamiltonian chaos using the Riemannian geometrization 
of newtonian dynamics. This work renewed an old intuition that dates back to N. S. Krylov 
11 1 and that spawned new ideas in abstract ergodic theory [|l2|,|l(|], whereas it did not give 



rise to any useful method to describe chaos in physical geodesic flows, despite of many 
attempts and with remarkable exceptions The obstacle was always the same: in 

analogy with Anosov flows, that live on hyperbolic manifolds, chaos has been invariably 
thought of as a consequence only of negative scalar curvature. So the first obvious check 
against any typical model that undergoes a stochastic transition - say the Henon-Heiles 
model - gives a puzzling surprise: the scalar curvature of (M,gj) is always positive || 
independently of the energy value, i.e. of regular or chaotic behavior of the dynamics. 

The novelty of the approach started in Ref. |J was to conjugate theoretical arguments 
with numerical experiments in order to shine some light on the following two points: i) does 
the geometry of the "mechanical" manifolds contain, though in some hidden way, the relevant 
information concerning stability and instability of their geodesies? and in the affirmative 
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case ii) how to quantify the strenght of chaos, how to characterize the weakly and strongly 
chaotic regimes? 

Actually positive answers to these questions have been given in []5| [TOj] , where, among 
other things, it has been shown that if the geodesies feel a positive non-constant curvature 
of the underlying manifold then parametric instability can be activated. Though a rigorous 
proof is not yet at disposal, parametric instability appears as the source of chaos on manifolds 
of positive non-constant curvature. 

By the way we can mention that also in the case of integrable systems, whose geodesies 
are therefore stable, the curvature of the underlying manifold can be wildly fluctuating along 
the geodesies but in this case the parametric instability mechanism is inactive, and it is found 
that these integrable geodesic flows have very special hidden symmetries, mathematically 



defined through Killing tensor fields ||32|| , that make them peculiar. 

For geodesic flows on constant negative curvature manifolds, the instability exponent 
is known [Eq. (pop], if the curvature is negative and non-constant then simple averaging 
algorithms can be devised, but what can we do with a positive and fluctuating curvature? 
The challenge was now to compute the average instability exponent for geodesic flows of 
physical relevance. This is a crucial test of effectiveness of the Riemannian theory of chaos 
with respect to the conventional explanation based on homoclinic intersections. Moreover, 
as no analytic method was available to compute Lyapunov exponents, it was worth making 
an effort in this direction. 

Under reasonable hypotheses, that obviously restrict the domain of validity of the ana- 
lytic formula ([)]]) for Ai, this paper provides the first analytic computations of the largest 
Lyapunov exponent in dynamical systems described by ordinary differential equations. 

Though several points need a deeper understanding, we hope that our work convincingly 
shows that this geometric approach is effective and useful, thus deserving further improve- 
ments and developments. 
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APPENDIX A: SOLUTION OF THE STOCHASTIC OSCILLATOR EQUATION 

In the following we will briefly describe how to cope with the stochastic oscillator problem 
which we encountered in Sec. [1 1 The discussion follows closely Van Kampen (Ref. ) 
where all the details can be found. 

A stochastic differential equation can be put in the general form 

F(x,Q) = 0, (Al) 

where F is an assigned function and the variable Q is a random process, defined by a mean, 
a standard deviation and an autocorrelation function. A function £ (Q) is a solution of this 
equation if VT2 F(£(Q), Q) = 0. If equation ( |A1|) is linear of order n, it is written as 

ii = A(t,Q)u (A2) 

where u6K" and A is a n x n matrix whose elements are randomly dependent on time. 

For the purposes of our work we are interested in studying the evolution of the average 
carried over all the realizations of the process , (u(t)). Let us consider the matrix A as the 
sum 

A(t,Q) = A (t) + aA 1 (t,tt) (A3) 

where the first term is ^-independent and the second one is randomly fluctuating with zero 
mean. Let us also assume that A is time-independent. If the parameter a - that determines 
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the fluctuation amplitude - is small we can treat Eq. ([AJ^) by means of a perturbation 
expansion. It is convenient to use the interaction picture, thus we put 



u(t) = exp(A t)v(t) 
Ai(i) = exp(A t)v(t) exp(— A t). 



(A4) 
(A5) 



Formally one is led to a Dyson expansion for the solution v(t). Then, going back to the 
previous variables and averaging, the second order approximation gives 
d f +co 

-{u(t)) = {A + a 2 j (A 1 (t)exp(A r)A 1 (t-r))exp(-A r)rfr}(u(t)) . (A6) 

Following the same procedure one can find also the evolution of the second moments (and 
by iterating also the evolution of higher moments). In fact, with the components of u G M n 
we can make n 2 quantities u v u^ that obey the differential equation 

d , 



(A7) 



where 



The above presented averaging method can be now applied to this new equation. 
Now, if we consider a random harmonic oscillator, Eq. (IA~3) has the form 



(A8) 
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with the random squared frequency Q = Qq + a n t](t). In particular, we are interested in 
working out the second moments equation when the process rj(t) is gaussian and 5-correlated. 
Using Eq. (|A8|) one finds that 



dt 





( 








2 


\ 


fxA 




I* 2 ) 


x 2 










-2Q 




x 2 


= A 


x 2 


, XX j 


\ 


-Q 


1 





/ 


, xx i 




, xx i 



(A10) 
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Because of our assumptions for this system, Eq. (|A6|) is more than a second order approxi- 
mation, it is exact. In fact, the Dyson series can be written in compact form as 



(x 2 (t)} 
\ (x{t)x{t)) J 



Kexp / A(t')dt' )1 



/ (^ 2 (o)) ^ 

<* 2 (0)> 
^ (x(O)x(O)) J 



(All) 



where the brackets [. . .] stand for a chronological product. According to Wick's procedure 
we can rewrite Eq. (|A11| ) as a cumulant expansion, and when the cumulants of order higher 
than the second vanish (as is the case of interest to us) one can easily show that Eq. (|A6|) 
is exact. 



Likewise in Eq. ( |A3| ), the matrix A splits as 



A(t) = A + a n 77(t)A 1 
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(A12) 



therefore the equation for the averages becomes 



d_ 
~dt 



1 (x 2 ) N 
(x 2 ) 

^ ( xx ) J 



+oo 



{A + a 2 / ( V (t) V (t - r))B(r)dr} 



1 (x 2 ) ^ 

(A 2 ) 
^ (xx) J 



(A13) 



where B(r) = Ai exp(A r)Ai exp(— A r). As (r)(t)r)(t — T)) = tS(t), with r a characteristic 
time scale of the process, we obtain 



' (x 2 ) \ ( (x 2 ) ^ 



d 
It 



(x 2 ) 
\ (xx) J 



{A + a 2 rB(0)} 



(x 2 ) 
\ (xx) J 



(A14) 



From the definition of B(r) it follows that B(0) = A 2 , then by easy calculations we find 
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(A15) 



which is the result used in Sec. 1IIB. 
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FIGURES 

FIG. 1. Average Ricci curvature (kn) vs. energy density e for the FPU model: comparison 
between analytic computation with Eq. ( [T4j ) (solid line) and the outcome of numerical simulations 
(time averages) with N = 128 (solid circles) and N = 512 (solid triangles); fi = 0.1. 

FIG. 2. Fluctuation of Ricci curvature (5 2 Kr)/N vs. energy density e for the FPU model: 
comparison between analytic computation with Eq. ( p2| ) (solid line) and numerical results. Symbols 
and parameters as in Fig. |l|. 

FIG. 3. Lyapunov exponent Ai vs. energy density e for the FPU model: comparison between 
theoretical prediction of Eq. ( |5l| ) (solid line) and numerical estimates at N = 256 (solid circles) 
and N = 2000 (solid squares); fi = 0.1. 

FIG. 4. Average Ricci curvature (kji) vs. energy density e for the coupled rotators model: 
comparison between analytic computation with Eq. (|92|) (solid line) and the outcome of numerical 
simulations (time averages) with N = 150 (solid circles); J = 1. 

FIG. 5. Fluctuation of Ricci curvature (5 2 Kr)/N vs. energy density e for the FPU model: 
comparison between analytic computation with Eq. (|97| ) (solid line) and numerical results. Symbols 
and parameters as in Fig. 

FIG. 6. Lyapunov exponent Ai vs. energy density e for the coupled rotators model: comparison 
between theoretical prediction of Eq. ( |5i~D (solid line) and numerical estimates at iV = 150 (solid 
circles), iV = 1000 (solid rhombs) and iV = 1500 (solid square); J = 1. 

FIG. 7. Estimate of the probability P (e) of occurrence of negative sectional curvatures in the 
coupled rotators model according to Eq. ( |101|) ; J = 1. 
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FIG. 8. Lyapunov exponent Ai vs. energy density e for the coupled rotators model: comparison 
between theoretical prediction and numerical estimates as in Fig. ||, but here the average curvature 
(fcft) which enters Eq. (|5l| ) is corrected according to Eq. (|102| ) with a = 150. Numerical values 
of Ai are obtained at N = 150 (solid circles), at N = 1000 (solid rhombs) and at N = 1500 (solid 
square) . 
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FIG. 2 Casetti et al. 

"Riemannian theory of Hamiltonian chaos..." 
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FIG. 3 Casetti et al. 
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FIG. 4 Casetti et al. 
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FIG. 5 Casetti et al. 
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FIG 6 Casetti et al. 
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FIG. 7 Casetti et al. 
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FIG. 8 Casetti et al. 
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